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Abstract.  We  present  a  reduced  basis  method  for  parametrized  partial  differential  equations  certified 
by  a  dnal-norm  bonnd  of  the  residnal  computed  not  in  the  typical  finite-element  “truth”  space  but 
rather  in  an  infinite-dimensional  function  space.  The  bound  builds  on  a  finite  element  method  and 
an  associated  reduced-basis  approximation  derived  from  a  minimum-residual  mixed  formulation.  The 
offline  stage  combines  a  spatial  mesh  adaptation  for  finite  element  and  greedy  parameter  sampling 
strategy  for  reduced  basis  to  yield  a  reliable  online  system  in  an  efficient  manner;  the  online  stage 
provides  the  solution  and  the  associated  dual-norm  bound  of  the  residual  for  any  parameter  value  in 
complexity  independent  of  the  finite  element  resolution.  We  assess  the  effectiveness  of  the  approach 
for  a  parametrized  reaction-diffusion  equation  and  a  parametrized  advection-diffusion  equation  with  a 
corner  singularity;  not  only  does  the  residual  bound  provide  reliable  certificates  for  the  solutions,  the 
associated  mesh  adaptivity  significantly  reduces  the  offline  computational  cost  for  the  reduced-basis 
generation  and  the  greedy  parameter  sampling  ensures  quasi-optimal  online  complexity. 
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1.  Introduction 

The  certified  reduced  basis  method  has  emerged  as  a  promising  model-reduction  strategy  to  rapidly  and 
reliably  solve  parametrized  partial  differential  equations  (PDEs)  in  real-time  and  many-query  scenarios  [13]. 
The  key  to  the  reliability  of  the  reduced-basis  method  is  a  posteriori  error  bounds.  However,  to  our  knowledge, 
with  the  few  exceptions  to  be  noted  shortly,  the  existing  bounds  are  with  respect  to  some  finite-element  “truth,” 
which  is  assumed  to  be  sufficiently  accurate  and  is  assumed  to  be  representative  of  the  infinite-dimensional 
weak  solution  of  the  PDE.  The  assumption  is  not  rigorously  validated  and,  especially  for  problems  with  spatial 
singularities,  may  be  violated  in  practice.  The  lack  of  reliable  a  posteriori  error  bounds  on  the  “truth”  leads 
to  either  an  unreliable  reduced-basis  approximation  with  respect  to  the  true  infinite-dimensional  weak  solution 
or  unnecessarily  expensive  computation  of  very  rich  finite-element  snapshots  in  the  offline  stage.  In  this  work, 
we  develop  a  reduced-basis  method  which  builds  a  residual  bound  relative  to  the  true  infinite-dimensional  weak 
solution;  we  hence  aim  to  remove  the  issue  of  truth  within  the  reduced-basis  framework. 
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The  residual  bound  strategy  proposed  in  this  work  is  based  on  a  minimum-residual  mixed  formulation  [4,10,11] 
of  parametrized  elliptic  PDEs.  Our  method  provides,  in  the  online  stage,  uniform  (as  opposed  to  asymptotic) 
bounds  for  the  dual  norm  of  the  residual  computed  relative  to  the  infinite-dimensional  function  space  for 
any  parameter  value.  We  emphasize  that  the  online  bound  is  1)  compnted  in  complexity  independent  of  the 
underlying  mixed  finite-element  discretization,  and  2)  available  for  any  parameter  value,  not  just  values  which 
belong  to  the  training  set  used  in  the  offline  stage.  The  proposed  minimum-residual  mixed  formulation  provides 
an  approximation  of  the  primal  held  as  well  as  the  dual  held,  the  latter  of  which  then  serves  to  form  a  bound 
of  the  dual  norm  of  the  residual.  In  addition,  like  other  minimum-residual  mixed  methods  and  unlike  Galerkin 
mixed  methods,  the  formulation  leads  to  coercive  hnite-element  and  reduced-basis  systems,  which  circumvent 
stability  issues  associated  with  for  instance  the  discrete  Ladyzenskaja-Babuska-Brezzi  (LBB)  condition;  the 
guaranteed  stability  allows  us  to  construct  the  hnite-element  and  reduced-basis  spaces  solely  from  approximation 
considerations. 

We  note  that  the  use  of  the  dual  variable  for  a  posteriori  error  estimates  —  in  particular  to  invoke  the  com¬ 
plementary  variational  principle  —  has  been  explored  in  the  past.  Some  of  the  earlier  works  in  the  hnite-element 
community  include  Ladeveze  and  Leguillion  [8],  Ainsworth  and  Oden  [1],  and  Sauer-Budge  et  al.  [15].  More 
recent  works  in  the  model-reduction  community  include  the  application  to  the  proper  generalized  decomposition 
by  Ledeveze  and  Chamoin  [7]  and  our  application  to  the  reduced-basis  method  [19].  In  particular,  in  terms  of 
the  objective,  our  earlier  work  [19]  is  similar  to  the  current  work  in  that  they  both  seek  a  strict  offline-online 
decomposition  of  the  exact  bound  computation  for  parametrized  PDEs.  However,  in  terms  of  the  formulation, 
our  earlier  work  [19]  differs  from  the  current  work  in  that  it  is  based  on  the  reduced-basis  approximation  of 
the  dual  held  that  satishes  exactly  the  so-called  dual  feasibility  condition  of  the  complementary  variational 
principle.  The  formulation  of  our  earlier  work  [19]  suffers  from  two  dehciencies:  hrst,  the  exact  satisfaction  of 
the  dual  feasibility  condition  requires,  in  the  offline  stage,  a  non-standard  hnite-element  approximation  and,  in 
the  online  stage,  a  potentially  large  algebraic  expansion  when  the  number  of  affine  expansion  terms  is  large; 
second,  the  application  is  limited  to  coercive  PDEs. 

In  order  to  overcome  the  limitations  of  the  bound  formulations  based  on  the  complementary  variational 
principle,  we  seek  in  this  work  the  dual-norm  bound  of  the  residual  relative  to  the  inhnite-dimensional  function 
space.  We  recall  that  the  dual  norm  of  the  residual  is  well-dehned  even  for  non-coercive  PDEs  and  is  equivalent 
to  the  norm  of  the  error.  Steih  and  Urban  [17]  have  recently  considered  the  problem  in  the  context  of  the 
reduced  basis  method;  their  method,  hrst,  utilizes  the  adaptive  wavelet  method  to  compute  the  reduced-basis 
snapshots  that  meet  the  user-specihed  true  error  tolerance  and,  second,  invokes  adaptive  residual  evaluation  in 
the  greedy  procedure.  We  note,  however,  the  adaptive  online  evaluation  of  the  dual  norm  of  the  residual  for  an 
arbitrary  parameter  value  requires  access  to  the  underlying  adaptive  wavelet  method  and  hence  is  not  online- 
efficient;  the  cost  of  the  online  certihcation  depends  on  the  complexity  of  the  underlying  wavelet  discretization. 
We  reiterate  that  our  goal  is  to  provide  an  online  bound  that  is  1)  computed  in  the  complexity  independent 
of  the  underlying  mixed  hnite-element  discretization,  and  2)  available  for  any  parameter  value,  not  just  values 
which  belong  to  the  offline  training  set. 

The  contributions  of  this  work  are  hvefold.  First,  we  develop  a  mixed  formulation  for  parametrized  elliptic 
equations  whose  dual  variable  naturally  yields  an  upper  bound  of  the  dual  norm  of  the  residual  associated 
with  the  primal  variable.  Second,  we  introduce  an  associated  mixed  hnite-element  and  mixed  reduced-basis 
approximations  for  the  weak  statement;  the  approximation  is  based  on  the  standard  Lagrange  and  Raviart- 
Thomas  elements,  which  facilitates  the  implementation  within  the  existing  hnite-element  libraries.  Third,  we 
develop  an  offline-online  computational  decomposition  of  the  mixed  formulation  for  parametrized  PDEs;  we  pay 
particular  attention  to  the  offline  construction  of  the  least-squares  mixed  system  in  the  parametrized  setting. 
Fourth,  we  introduce  a  spatio-parameter  adaptation  strategy  based  on  an  isotropic-/i  mesh  adaptation  in  space 
and  a  greedy  sampling  in  parameter.  Fifth,  and  hnally,  we  demonstrate  the  effectiveness  of  the  proposed  strategy 
applied  to  an  elliptic  problem  with  spatial  singularities. 

We  note  an  important  limitation  of  the  current  work.  While  we  provide  a  dual-norm  bound  of  the  residual 
which  is  equivalent  to  the  norm  of  the  error,  we  do  not  provide  a  lower  bound  of  the  stability  constant  which 
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is  required  to  quantitatively  bound  the  norm  of  the  error.  The  computation  of  a  lower  bound  of  the  stability 
constant  relative  to  the  infinite-dimensional  function  space  is  beyond  the  scope  of  this  work.  We  however  note 
that  residual  bounds  or  even  residual  estimates  have  often  been  found  to  be  sufficient  for  the  purpose  of  both 
spatial  mesh  adaptation  for  finite  element  [1]  and  greedy  parameter  sampling  for  reduced  basis  [5,18]. 

This  paper  is  organized  as  follows.  In  Section  2,  we  devise  a  residual-bound  strategy  based  on  the  approxi¬ 
mation  of  the  primal  and  dual  fields  and  introduce  the  associated  minimum-residual  mixed  finite-element  and 
reduced-basis  methods.  In  Section  3,  we  develop  an  offline-online  computational  strategy  for  the  minimum- 
residual  mixed  reduced-basis  method.  In  Section  4,  we  develop  a  spatial  mesh  adaptation  strategy  and  a 
greedy  parameter  sampling  strategy.  In  Section  5,  we  apply  the  method  to  a  reaction-diffusion  equation  and  an 
advection-diffusion  equation. 


2.  Bound  Construction 

2.1.  Problem  Description 

We  first  introduce  a  d-dimensional  bounded  domain  O  C  with  a  Lipschitz  boundary  90  such  that  Tat  U 
r £)  =  90  and  T^)  is  nonempty.  We  also  introduce  a  P-dimensional  bounded  parameter  domain  T)  C  We 
then  consider  parametrized  boundary  problems  of  the  following  form:  given  a  parameter  /i  €  P,  find  u(/r)  that 
satisfies 


V  •  Vu(^);  fj.)  +  C{u{n);  fi)  =  /(/r)  in  O, 

■n  =  on  Tat, 
u(/i)  =  0  on  Td; 

here  A^WjVw;  fj,)  is  a  flux  operator  that  is  linear  in  each  of  w  and  Vru,  C{w;g)  is  a  reaction  operator  that  is 
linear  in  w,  f{g)  is  a  source  function,  and  g{g)  is  a  Neumann  boundary  source  function.  We  will  shortly  state 
more  precise  conditions  on  the  flux,  reaction,  and  source  operators. 

In  order  to  cast  the  problem  in  a  weak  form,  we  now  introduce  a  Hilbert  space 

V  =  {v&H\^)  I  u|r„  =0} 

endowed  with  an  inner  product  {w,  v)v  =  Vw  ■  Vvdx  +  wvdx  +  wvds  and  the  induced  norm  ||rc||  v  = 
a/ {w,  rc)v-  We  then  consider  the  following  weak  problem:  given  g  gV,  find  u{g)  €  V  such  that 

a{u{g),v;  g,)  =  £{v,  fi)  Vu  G  V, 

where 


a{w,v;g)  =  —  /  S/v  •  A{w,'\/w;  g)dx  +  /  vC{w,g)dx, 
Jn  Jn 

£{v;  g)  =  f  vf{g)dx—  [  vg{g)ds  Vu  G  V. 


Jn  JFk 

We  require  that  H  :  V  x  {L‘^{Q)Y  x  I?  — >  (L^(H))'^,  C  :  V  x  I?  — >  /:!?—>  and  g  :  V  ^  L‘^{T at); 

we  assume  A{-,-;g),  C{-;g),  f{g),  and  g{g)  are  bounded  for  all  g  GV.  We  also  assume,  to  ensure  the  problem 
is  well-posed,  that  the  bilinear  form  a(-,  •;  g)  is  inf-sup  stable  and  continuous: 


o/  \  •  r  a{w,v;g) 

P{g)  =  mf  sup  t; — n — ,,  ,,  >  0 


'y{g)  =  sup  sup 


a{w,  v;  g) 


<  oo 


mevwev  ll^llvll^'llv 


VgGV, 
Wg  G  V. 


(1) 


(2) 
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In  addition,  we  assume  that  A,  C,  f,  and  g  permit  an  affine  decomposition  in  the  sense  that 


Qa 

A{w,Vw;fi)  ='^Qgi^i)Aq{w,Vw), 

9=1 

Qc 

C{w;  m)  =  X!  ®9(m)C'9(i«)> 

9=1 

Qf 

/(m)  = 

9=1 

Qg 

gip)  =  '^^q{g)9q, 

9=1 

for  some  parameter-independent  bounded  operators  Aq  :  Vx  :  V  — >  fq  € 

and  e  L^(r jv)  and  parameter-dependent  functions  0“  :  I?  — >  M,  0°  :  I?  — )■  M,  0^  :  I?  — )■  M,  and  0®  :  I?  — >•  M. 
We  note  that  this  is  the  so  called  “affine  decomposition”  assumption  that  is  common  to  many  reduced  basis 
methods  [13];  for  operators  that  do  not  provide  an  affine  decomposition,  the  empirical  interpolation  method  [2] 
may  provide  an  approximate  decomposition. 

We  now  define  the  residual.  For  any  w  €  V,  the  residual  evaluated  relative  to  some  u  €  V  is  given  by 


r{v;  w;  g)  =  £{v;  g)  —  a{w,  v;  g) 

=  /  vf{g)dx—  /  vg{g)ds  +  /  'S/v  ■  F{w^'S/w]  g)dx  — 

J  Vi  r  jv  J  Vi 

We  wish  to  construct  an  upper  bound  of  the  dual  norm  of  the  residual 

\r{v,w]g)\ 


vC{w;  g)dx. 


\'>'i-;w;g)\\v'  =  sup  - 


uev  Ipllv 


where  the  element  w  €  V  arises  from  a  reduced-basis  approximation.  We  note  that  the  residual  is  related  to 
the  error  \\u{g)  —  m||vj  for  any  w  €  V,  by 

||m(m)  -  w||v  <  ■g^||r(-;w;/9)||v'  <  ||w(ai)  -  w|| v, 

where  I3{g)  is  the  inf-sup  (or  stability)  constant  (1)  and  7(/r)  is  the  continuity  constant  (2).  In  this  sense, 
the  dual  norm  of  the  residual  is  equivalent  to  the  norm  of  the  error.  However,  in  practice,  a  lower  bound  of 
the  inf-sup  constant  is  required  to  quantitatively  bound  the  norm  of  the  error.  As  noted  in  the  Introduction, 
in  this  work  we  focus  on  the  computation  of  an  upper  bound  of  the  dual  norm  of  the  residual  and  devise  a 
computational  strategy  that  admits  an  offline-online  decomposition;  the  computation  of  a  lower  bound  of  the 
inf-sup  constant  is  beyond  the  scope  of  this  work. 

Remark  1.  Many  elliptic  equations  conform  to  the  form  introduced  above.  For  instance,  a  diffusion  equation  — 
with  a  parametrized  diffusion  coefficient  K{g)  €  M>o  —  results  from  A{w,Vw;  g)  =  —K{g)'Vw.  An  advection- 
reaction-diffusion  equation  —  with  parametrized  diffusion  coefficient  K{g)  €  M>o,  advection  coefficient  f3{g)  € 
and  reaction  coefficient  ^{g)  €  IR  —  results  from  A{w,  Vw;  g)  =  —K{g)S/w  +  (d{g)w  and  C{w,  g)  =  j{g)w. 
The  Helmholtz  equation  —  with  the  wavenumber  as  the  parameter  —  results  from  A(w,  Vm;  g)  =  — Vw  and 
C{w;  g)  =  —g^w.  The  form  also  supports  the  variants  of  the  problems  with  parametrized  fields  instead  of  the 
constants. 
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2.2.  Primal-Dual  Residual  Bouud:  Theory 

We  introduce  a  Hilbert  space 

Q  =  R(div;  O)  =  {<7  e  |  V  •  r  G  L^{n)} 

endowed  with  an  inner  product  (p,  q)Q  =  {p,  q)H{div,n)  =  JqP  '  +  /n(^  '  P)i'^  '  Q)dx  and  the  induced  norm 

llgljg  =  ||(7||Lr(div;a)  =  \/{q,  q)H(div;n)-  We  are  now  ready  to  state  the  key  proposition  of  this  work. 

Proposition  1.  For  any  p,  G  T>  and  w  G  V,  the  dual  norm  of  the  residual  is  bounded  by 

\\r{-;w;p)\\v'  <  ^/B{w,q;p)  Vq  G  Q, 

where 

B{w,q;p)  =  \\f{p)  -  C{w;p)  -  V  •  gf  -t  \\A{w,Ww;p)  -  gf  -k  \\gip)  -  9  •  n|li2(r„)-  (3) 

Proof.  For  notational  simplicity,  in  this  proof,  we  suppress  the  explicit  appearance  of  the  parameter  p  in  various 
forms  and  operators.  We  first  recall  that  the  dual  norm  may  be  expressed  in  terms  of  its  Riesz  representation: 

l!R(w)||v  =  lk(-;u;)||v'  Vw  G  V, 

where  the  Riesz  representation  R{w)  G  V  satisfies 

{R{w),v)v  =  t{v;w)  Vp  G  V. 

We  then  consider  the  following  inequality:  for  any  v  G  V, 


Vv  ■VR{w)dx  +  /  vR{w)dx  +  /  vR{s)ds 
I  Jn  Jfk 

/  vfdx—  /  vgds  +  /  S/v  ■  A{w,'\/w)dx  — 
«/  o  r  TV  *^0 

/  vfdx—  /  vgds  +  /  X^v  •  A{w,\/w)dx  — 

J  Q,  ^  r  AT  J  Cl 


—  /  vV  •  qdx  +  /  vq  '  fids  —  /  Vv  •  qdx 


I 

Jsi 

j 

JQ. 


vC{w)dx 

vC{w)dx 


JQ. 


f  v{f  —  C{w)  —  S/ ■  q)dx 


a 

v{g  —  q  ■  n)dx  +  (  \/v  •  {A{w,'\/w)  —  q)dx 


J  Cl  ^  r  AT  ^  Cl 

<  ll^'l|L2(n)ll/  -  C{w)  -  V  •  qU^Q)  +  l|v||L2(r„)||5  -  <? '  n\\mr^)  +  |!Vp||L2(n)||H(u;,  Vw)  -  q\\L^(Q) 

^  (lklli2(a)  +  ||v||i2(r„)  +  l|Vn|||2(Q)^ 

(ll/  -  C'M  -  V  •  qWh^a)  +  \\9-q-  «lli2(r«)  +  ^(w,  Vw)  -  qWh^ft))  ^ 

=  \\v\\viBiw,q))^/^; 


here  the  second  equality  follows  from  the  Green’s  theorem,  the  first  inequality  follows  from  the  triangle  inequal¬ 
ity,  the  second  equality  follows  from  Cauchy-Schwarz,  and  the  last  eqnality  follows  from  ||i;||y  =  ||p|li2(n)  + 
ll^llL2(rN)  II^^IIl2(o)  and  the  definition  of  the  form  B{-,  •).  We  finally  set  v  =  R{w)  and  invoke  the  equivalence 
||i?(n;)||v  =  ||r(-; rc)||v'  to  obtain  the  desired  result.  □ 
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Remark  2.  We  may  in  principle  consider  more  general  source  functions  —  namely  /(/i)  €  H  and  g(n)  G 
^-i/2(r^)  —  modifying  the  first  inequality  in  the  proof  to 

•  •  •  <  MmmWf  -  C{w)  -  V  •  g|lff-i(a)  +  ||v||Hi/2(r„)||5  -  q  ■  ^lli/-i/2(r„) 

+  l!Vp||L2(Q)||yl(w;,  Vw)  -  gWr^ii)- 

However,  the  counterpart  of  the  form  associated  with  this  decomposition  requires  the  and 

norms,  which  are  not  readily  computable.  We  hence  refrain  from  the  consideration  of  these  more 
general  source  functions  in  this  work. 

2.3.  Minimum-Residual  Mixed  Finite  Element  Method 

Our  reduced-basis  approximation,  as  in  the  standard  reduced-basis  approach  [13],  builds  on  the  finite-element 
snapshots  computed  at  select  parameter  values.  We  hence  first  present  the  finite-element  discretization  used 
in  this  work.  We  note  that  the  discretization  is  a  generalization  of  least-squares  mixed  finite  element  methods 
reviewed  by  Bochev  and  Gunzburger  [4]. 

We  first  introduce  a  triangulation  Th  oi  Q.  We  then  introduce  conforming  finite-element  approximation 
spaces  for  V  and  Q: 


V^^  =  {vGV\  vU  gW,  kg  %}, 

Q^^  =  {qGQ\  qU  G  ©  x¥P-\  k  G  %}; 

note  that  the  space  consists  of  Raviart-Thomas  elements  [12]  of  degree  p  —  1.  The  superscripts  Afv  and 
A/q  signify  the  number  of  degrees  of  freedom  associated  with  the  spaces  and  ,  respectively. 

We  now  consider  the  following  minimum-residual  approximation  problem:  given  p  G  V,  find  {u^ {p) ,  p^ {ii))  G 
X  Q-^®  such  that 


{u-^{p),p-^{p))  =  aiginf  B{w,q;p),  (4) 

wGV^v 

qGQ^a 

where  B{-,-;-)  is  the  bound  form  (3).  We  may  identify  the  Euler-Lagrange  equation  associated  with  the  mini¬ 
mization  statement:  given  p  GV,  find  {u^ {p),p^ {p))  G  x  Q-^®  such  that 

G{{u-^{p),p-^{p)),{v,q);p)  =  L{{v,q);p)  Vv  G  ,  Vq  G  ,  (5) 


where 


G{{w,r),{v,q);p)  =  [  {C{v;  p) +  V  ■  q){C(w;  p) +V  ■  r)dx 
Jn 

+  /  {A{v V,  p)  —  q)  ■  {A{w W]  p)  —  r)dx 

Jn 

+  /  {q  •  h){r  •  h)ds, 

Jfn 

L{{v,q)]p)=  /  {C{v,p)  +  \I  ■q)f{p)dx+  /  {q  ■  n)g{p)ds; 
Jn  Jtw 


(6) 


we  take  Af  =  Afv  +  A/g,  which  serves  as  a  measure  of  the  complexity  of  the  finite-element  approximation. 
We  have  the  following  proposition  as  regard  the  coercivity  of  the  bilinear  form  G(-,  •;  •). 
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Proposition  2.  The  bilinear  form  G(-,  •;  jj,)  defined  in  (6)  is  coercive  for  any  jj,  G  T>: 


+  Ugllg)  Vi;  G  V,  Vg  e  Q, 


for  a  coercivity  constant 


a{fj.)  = 


Wp)y‘ 


(l  +  2||G(-;  M)ll£(v,L2(n))  +  2||^(-,  •;  M)ll£(v,(L2(n))rf; 


here  is  the  inf-sup  constant  (1). 

Proof.  The  proof  is  provided  in  Appendix  A.  □ 

Because  the  bilinear  form  is  symmetric,  bounded  (by  inspection),  and  coercive,  the  problem  (5)  is  well-posed. 
Note  in  particular  the  minimum-residual  mixed  finite  element  method  is  not  subject  to  the  LBB  condition,  as 
noted  earlier  for  coercive  PDEs  by  Pehlivanov  et  al.  [10, 11]  and  discussed  in  the  review  for  other  classes  of 
PDEs  by  Bochev  and  Gunzburger  [4]. 

The  residual  associated  with  our  mixed  finite-element  approximation  u^{p,)  G  V'^'^  is  bounded  by,  for  any 
fj,  GP, 


In  addition,  because  the  bound  functional  coincides  with  the  objective  function  of  the  minimum-residual  state¬ 
ment,  the  bound  decreases  monotonically  in  the  sense  that 

B{u^'  <  B{u^{p,),p-^{p,);p) 

for  nested  finite-element  spaces  D  and  D  We  will  later  appeal  to  this  monotonicity  to 

develop  a  convergent  adaptive  algorithm. 

2.4.  Minimum-Residual  Mixed  Reduced  Basis  Method 

We  now  introduce  our  reduced-basis  approximation.  We  hrst  introduce  an  A^-dimensional  primal  approxi¬ 
mation  space  spanned  by  a  basis  G 

Vat  =  span{^„}()Li  C 

and  an  iV-dimensional  dual  approximation  space  spanned  by  a  basis  {?7„  G 

Qn  =  span{?7„}()Li  C  Q"^®. 

We  then  identify  the  reduced-basis  counterpart  of  the  minimum-residual  hnite-element  statement  (4):  given 
fj,  GT),  hnd  (MAr(/x),pAr(/x))  G  Vn  x  Qn  such  that 

{uN{p),PNi^J'))  =  argM  B{w,q;p).  (7) 

wGVjv 

The  associated  reduced-basis  Euler-Lagrange  equation  —  the  counterpart  to  the  finite-element  statement  (5) 
—  is  as  follows:  given  p  GV,  hnd  {un{p),Pn{p))  G  Vn  x  Qm  such  that 


G{{un{p),Pn{p)),  (v,  q);M)  =  L((v,  q);p)  Vw  G  Viv,  Vg  G  Qn, 


(8) 
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where  the  bilinear  form  and  linear  form  are  as  defined  in  (6).  Because  the  bilinear  form  is 

symmetric,  bonnded,  and  coercive,  the  reduced-basis  system  is  well-posed.  Note  in  particular  the  minimum- 
residual  mixed  reduced  basis  method,  as  in  the  finite-element  counterpart,  is  not  subject  to  the  LBB  condition; 
this  is  unlike  many  Galerkin  mixed  reduced  basis  methods,  which  require  a  careful  choice  of  the  spaces  (c.f.  [14]). 

Analogous  to  the  finite-element  case,  the  residual  associated  with  our  mixed  reduced-basis  approximation 
un{^)  €  Vat  is  bounded  by 


In  addition,  because  the  bound  functional  coincides  with  the  objective  function  of  the  minimum-residual  state¬ 
ment,  the  bound  decreases  monotonically  in  the  sense  that 

B{uN'{p),PN'{p);p)  <  b{un{p),pn{p)\p) 

for  hierarchical  reduced-basis  spaces  Vn'  3  Vn  and  Qn'  3  Qn-  We  will  later  appeal  to  this  monotonicity  to 
develop  a  convergent  adaptive  algorithm. 

Remark  3.  Because  the  bilinear  form  associated  with  the  minimum-residual  statement  is  coercive,  we  have 
great  flexibility  in  choosing  the  reduced-basis  approximation  spaces.  For  instance,  we  could  consider  primal 
and  dual  approximation  spaces  of  different  dimensions.  As  an  another  example,  instead  of  considering  separate 
reduced-basis  spaces  for  the  primal  and  dual  variables,  we  could  couple  the  primal  and  dual  basis  functions  and 
consider  a  reduce-basis  approximation  of  the  couple;  the  coupled  approach  in  general  results  in  a  lower  online 
cost  but  a  higher  offline  cost. 


3.  Computational  Considerations 

3.1.  Assumption:  Polynomial  Form 

We  here  develop  an  efficient  computational  strategy  for  a  special  case:  the  flux  operator  A(-,  •;  p)  and  reaction 
operator  C{-;  p)  preserve  the  piecewise  polynomial  structure  of  the  approximation  space  V'^^;  the  source  terms 
f{p)  and  g{p)  are  piecewise  polynomial  consistent  with  the  triangulation  Th-  More  precisely,  we  hrst  introduce 
element-wise  degree-A:  polynomial  discontinuous  approximation  spaces 

=  {xG  L^{n)  I  x\,  e  p^  Vk  e  %}, 

=  {ye  \  y\,  e  Vk  g  %}, 

Z^^={zGL\r)  I  ^UGP^  aGS4, 

for  'Eh  the  surface  triangulation  of  Fat  that  is  consistent  with  the  volume  triangulation  7h  of  P;  we  then  require 
that,  for  all  p  GV, 


A{w,  Ww;  p)  G  y^y  Vw  G  , 

C{w;p)gX^^ 

f{p)GX^y 

g{p)&Z^-. 

Note  that  the  polynomial  degree  k  need  not  be  the  same  as  the  polynomial  degree  of  the  space  p;  in 
particular  we  may  consider  k  >  p  to  treat  a  larger  class  of  operators. 

We  emphasize  that  the  finite-element  and  reduced-basis  formulations  presented  in  the  previous  section  applies 
to  general  A(-,-;/x),  C{-;p),  f{p),  and  g{p);  the  formulations  do  not  require  the  operators  to  preserve  the 
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piecewise  polynomial  structure.  We  may  also  readily  derive  an  associated  computational  strategy  for  the 
general  operators.  However,  when  the  operators  preserve  the  piecewise  polynomial  structure,  we  can  greatly 
simplify  the  implementation  of  the  minimum-residual  mixed  finite  element  and  reduced  basis  methods  and 
in  addition  reduce  their  computational  complexities.  We  also  note  that  many  flux  and  reaction  operators  of 
practical  interest  preserve  the  polynomial  structure. 

3.2.  Discrete  Operators 

We  now  introduce  discrete  operators  (i.e.  matrices)  that  will  be  used  in  our  finite-element  and  reduced- 
basis  approximations.  We  first  introduce  a  basis  {Xm\m=i  of  and  define  an  element-wise  block  matrix 

X  G  M-'^vxa/'a-  entries 

Xmn=  [  XmXndx. 

Ja 

We  next  introduce  a  basis  {7m}m=i  of  and  define  element-wise  block  matrix  Y  e  with  entries 

^mn  =  /  Tm  *  ^ndx. 

Jn 

We  then  introduce  a  basis  define  element-wise  block  matrix  Z  G 


^m.n,  — 


we  note  that  the  matrices  X,  Y,  and  Z  are  symmetric  positive  definite.  We  now  introduce  a  basis  of 

and  define  matrices  R  G  p)  ^  ]gY.Yx7VQ^  g  ^  entries 

R-mn  =  /  7m  ‘ 

JQ 

^mn  =  /  Xm{^  * 

JQ 

^mn  =  /  Cm‘4^n  *  uds. 

We  finally  introduce  a  basis  of  and  define  matrices  A(/i)  €  C(/i)  €  R-^vxYv  g^^j 

vectors  f(/i)  G  and  g(/i)  G  R-^^  with  entries 

Qa 

I 

9=1 


p  p  '‘Hi  a  a 

Am„(/x)  =  /  7m  •  A((/)„,  =  /  7m  •  ^  ^  e“(^)A^„, 

Cm„(/l)=  f  XmC{^n\P^)dx=  f  Xm  ^  0^(/r)C^„, 

^  „  Qf  Qf 

fm(/r)  =  /  Xmf{fJ^)dx=  /  = 

g=l  9=1 

^  p  Qa  Qa 

gm(/l)  =  /  Cm9ia')dx  =  Cm 
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where  A"?  €  ,  q  =  1, . . . ,  Qa,  ^  q  —  f?  g  ,  q  =  1,...,(5/,  and  g'^  g  , 

g  =  1, . . . ,  Qg,  are  given  by 


Ain=  /  7m  •  V(/)„)dx, 

Jq 

Cl„=  [  XmCH^n)dx, 

Jq 

fm  =  /  Xm^dx, 

Jq 

gm  =  /  Cmff'^t^s; 

Jr„ 

note  that  the  construction  appeals  to  the  affine  decomposition  of  the  operators. 


3.3.  Implementation  of  the  Minimum-Residual  Mixed  Finite  Element  Method 

We  hrst  note  that,  for  piecewise-polynomial  preserving  A(-,-;/i)  and  C'(-;/r),  the  residual  bound  operator 
maybe  expressed  as,  for  all  w  €  V'^'^  and  q  € 

B{w,  q;  ^l)  =  ||/(/r)  -  C{w]  n)-V  ■  +  ||  A(u;,  Vw;  ^l)  -  qWl^^a) 

+  Wdid-)  —  Q  ■  '^IIl2(p^) 

=  l|nA'V;^[/(/i)  -  C{w;ii)-\/  ■  g]||p2(n)  +  \\Ilyj~fy[A{w,S/w,  fi)  -  g]||i2(n) 

+  Iln^Yz  [^(/x)  —  q  ■  h]||p2(pj,^), 

where  ,  and  are  the  orthogonal  projection  operators  onto  ,  and  ,  respectively. 

We  may  thus  express  the  residual  bound  B{-,  y/i)  dehned  in  (3)  using  the  discrete  operators  as  follows:  for  any 
■l"  =  J2m=l  WTn</>m  and  q  =  J2mil  hmV'm, 


where 


B{w,q;fi) 


||f(/x)  -  C{fi)w  -  Dq|l|-_i  +  ||A(/x)w  -  RqUy-i  +  ||g(Ax)  -  Bq|l|-i 


(qV/KM  KM 

\  -W  J  y  K2i(/x)  K22(Ai) 


q  \  _2  Li(/x)  q 

w  y  L2(/x)  J  y  w 


+  F(/i) 


Kii(^)  =  D^X'^D  +  +  B^Z’^B, 

Ki2(ax)  =  K^i(ax)  =  D^X-iC(/x)  -  R^Y”iA(/x), 

K22(m)  =  +  A(/x)^Y-1A(ax), 

M)  =  D^X-if(Ax)  +  B^Z-ig(M), 

F{fi)  =  f(/x)^X-if(Ax)  +  g(M)^Z-ig(Ax). 

The  discrete  Euler-Lagrange  equation  is  the  following:  find  (p(/x)  €  K-^'^,u(/x)  G  K-^s)  such  that 

f  Kn(/x)  Ki2(m)  \  f  p(m)  Li(m)  \  . 

K2i(/x)  K22(m)  J  \  u(/i)  J  y  M)  J  ’ 
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In  practice,  we  first  compute  the  linear  combination  of  the  discrete  operators  to  form  A(/x)  = 

C(m)  =  f(/x)  =  E^=i  and  g{^i)  =  &gifJ.)gq-  We  then  form  the  matrices  (^), 

Li(/r),  and  F(^)  defined  above.  We  finally  solve  the  linear  system.  This  construction  of  the  discrete  system  — 
only  permitted  for  the  piecewise-polynomial  preserving  operators  —  requires  computation  of  0(max{Q°' ,  Q'^}) 
operators  and  in  particular  avoids  the  computation  of  the  0{max{Q°',Q‘^}'^)  operators  that  would  arise  from 
a  direct  computation  of  the  least-squares  operators.  Note  that  the  latter  can  be  prohibitive  for  systems  that 
require  a  large  number  of  affine  expansion  terms. 

Remark  4.  For  operators  that  that  do  not  preserve  the  piecewise-polynomial  structures,  we  may  obtain  the  dis¬ 
crete  residual  and  the  discrete  Euler-Lagrange  equations  through  a  direct  substitution  of  the  basis-representation 
of  functions  into  (3)  and  (5),  respectively. 


3.4.  Implementation  of  the  Minimum-Residual  Mixed  Reduced  Basis  Method 

We  now  consider  our  reduced-basis  implementation.  We  first  note  that  we  may  express  any  w  €  Vat  C 
as  w  =  rnn4‘mOLn  for  some  coefficient  vector  a  G  and  a  matrix  V  G  rNvxm 

the  n-th  reduced  basis  is  Similarly,  we  may  express  any  q  G  Qn  C  as  q  = 

EmSi  Qmn'ipmfin  fo''  some  Coefficient  vector  (3  G  and  a  matrix  Q  G  such  that  the  n-th  reduced 

basis  is  =  Em£i  Qmn7n-  We  may  thus  express  the  residual  bound  B(-,  •;  fj,)  specialized  to  the  reduced-basis 
spaces  using  the  discrete  operators:  for  any  w  =  J2n=i  Em=i  mn4'mOLn  and  q  =  J2n=i  EmSi  Qmni’mf^n, 


B{w,q;  fi) 


Y  (  Q^Kn(Ai)Q 
)  V^K2i(m)Q 

Q^Li(ai)  /3 
V^L2(/x)  )  \o, 


Q^Ki2(ai)V 

V^K22(/x)V 

)  +F(ai). 


/3 

a 


The  discrete  reduced-basis  Euler-Lagrange  equation  is  the  following:  hnd  such  that 


Q^Kn(M)Q  Q^Ki2(m)V  \  ( 
V^K2i(/x)Q  V^K22(m)V  J  [ 


Q^Li(ai) 

V^L2(/x) 


Here  we  provide  explicit  representations  of  the  reduced-basis  matrices  using  the  affine  expansion: 
Q^KiiQ  =  [Q^D^X-^DQ]  +  [Q^R^Y-^RQ]  +  [Q^B^Z-^BQ]  , 

Q^Ki2V  =  ^  0=(/x)  [Q^D^X-iC«V]  -  ^  0“(m)  [Q^R^Y-1A«V]  , 

9=1  9=1 

Qo  ^  ,  Qa 


V^K22V=  ^  0^,(m)0^(m)  [v^c^'x-^c^vj -f  ^  0“,(a^)0“(m) 

9', 9=1  9'. 9=1 

Q  f  Qg 

Q^Li(m)  =  £  0^(/i)  [Q^D^X-iF]  +  £  09(^)  [B^Z-ig«]  , 

9=1  9=1 

Qc  Qf 

^  ^0^(M)e^(Ai)  [(C«yx-1F 

9'  =  1  9=1 


yi  A«  Y-^A«V 


The  offline-online  computational  decomposition  is  clear  from  the  construction.  In  the  offline  stage,  we  com¬ 
pute  the  finite-element  coefficients  associated  with  the  reduced  basis,  V  and  Q,  and  then  form  the  various 
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parameter-independent  matrices  and  their  parameter-independent  products  indicated  by  brackets  in  the  ex¬ 
pression  above.  In  the  online  stage,  we  hrst  take  appropriate  linear  combinations  of  the  parameter-independent 
matrices  weighted  by  parameter-dependent  functions  to  form  the  reduced-basis  Euler-Lagrange  system  of  size 
2N  X  2N-,  we  then  solve  the  reduced-basis  system  for  the  coefficient  vectors  Q.*{y)  and  /3*(/r). 

4.  Spatio- Parameter  Adaptation 

4.1.  Adaptation  of  the  Finite  Element  Spaces 

Although  ultimate  goal  is  to  provide  the  reduced-basis  spaces  Vn  and  Qn,  approximations  in  which  meet  a 
user-specified  residual  tolerance,  we  first  briefly  discuss  our  hnite-element  adaptation  procedure  that  generates 
and  because  Vn  C  and  Qn  ^  by  construction.  In  this  work,  we  consider  for  simplicity 

isotropic  h  adaptivity  for  a  fixed  approximation  order  p;  however,  the  residual  bound  procedure  can  readily 
accommodate  general  anisotropic  hp  finite-element  approximation  spaces. 

We  introduce  an  elemental  residual  indicator  associated  with  our  hnite-element  approximations  {w^ {p,),p^{p))  € 
yNv  X  qAq. 


we  observe  that 


B{u^{p),p^{p);p)  =  rj^. 

Once  we  identify  largest  contributors  to  the  residual  bound,  we  employ  a  hxed-fraction  marking  strategy  to 
mark  elements  to  rehne;  in  our  particular  work  we  rehne  the  top  10%  of  the  elements  with  the  largest  residual 
indicator.  We  then  rehne  the  marked  elements  using  the  newest-vertex  bisection  algorithm  [9];  we  recall  that 
the  algorithm  creates  strictly  nested  meshes  such  that  C  C  . . .  and  C  C  . . .  for  a  sequence 
of  adapted  meshes.  The  construction,  together  with  the  fact  that  the  solution  pair  is  the  minimizer  of  the 
residual  bound,  implies  that  the  residual  bound  monotonically  decreases  with  each  mesh  rehnement. 

4.2.  Construction  of  Vat  and  Qn-  Spatio-Parameter  Greedy  Algorithm 

We  construct  the  reduced-basis  spaces  Vn  and  Qn  using  a  spatio-parameter  greedy  procedure  described  in 
Algorithm  1.  We  here  describe  the  algorithm  in  words.  We  first  identify  the  parameter  whose  state  is  least- 
well  represented  in  the  reduced-basis  spaces  Vat  and  Qat  as  assessed  by  the  residual  bound  with  respect  to 
the  infinite-dimensional  function  space.  We  then  compute  the  associated  finite-element  approximation  of  the 
primal  and  dual  fields  using  the  current  finite-element  discretization,  V^^  and  .  If  the  residual  certificate 
B{w^'‘{flN),P^’‘ipN)',  P^n)  associated  with  the  finite-element  solution  is  greater  than  the  prescribed  tolerance, 
then  we  adaptively  refine  the  hnite-element  space  until  we  meet  the  desired  tolerance;  recall  that  the  residual 
certihcate  decreases  monotonically  with  the  nested  rehnement  of  the  hnite-element  space,  hence  we  are  guar¬ 
anteed  to  meet  the  tolerance  upon  sufficient  rehnement.  We  then  augment  the  reduced-basis  spaces  with  these 
hnite-element  solutions  that  meet  the  required  error  tolerance. 

We  make  a  few  remarks.  First,  if  the  initial  hnite-element  approximation  spaces  are  sufficiently  rich  in 
the  sense  that  sup^gg^^^.^  inf^gyAfo  B{w,q;p)  <  then  the  algorithm  works  as  the  standard  WeakGreedy 

q&Q^O 

algorithm  [13]:  the  hnite-element  spaces  are  not  enriched  as  the  rehnement  is  not  necessary  to  meet  the  desired 
tolerance.  Second,  whenever  the  hnite-element  space  is  enriched,  we  re-represent,  in  terms  of  the  computational 
implementation,  the  reduced  basis  computed  on  an  earlier  mesh  on  the  most  recent  mesh;  mathematically,  the 
reduced  basis  are  unaltered  since  the  strictly  nested  meshes  imply  that  the  reduced  basis  computed  on  the 
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Algorithm  1:  Spatio-parameter  Greedy  algorithm 
input  ;  ,  Q^° :  initial  finite-element  spaces 

etoi  €  residual  tolerance 

output:  Swmax-  reduced-basis  approximation  spaces 

1  for  A  =  1,  ,  Amax  do 

2  Identify  the  parameter  associated  with  the  largest  residual  bound 

/iTV  =  argsup  inf  B{w,q;^); 
weVN 


if  inf.a,eVjv  B{w,  q;  p)  <  e^^i,  terminate. 

qG  Qn 

3  Evaluate  the  associated  finite-element  approximation 

=  inf  B{w,q]flN) 

qGQ^Ic 


4  if  then 

5  Enrich  the  finite-element  spaces  k'  times  to  meet  etoi  for  /In 

pAfc  ^  gAfc  pAfc+fc/  ^  y^^k  ^  QJ^k+y  ^  gAfc  ^ 


6  Set  k  k  +  k' . 

7  Augment  reduced-basis  spaces 


Viv  =  span{VAr_i,if^'=(/i7v)}  and  Qiv  =  spanlQjv-uP'^HAw)}- 


8  end 


earlier  mesh  is  exactly  representable  on  all  subsequent  meshes.  Third,  the  algorithm  is  convergent  since  the 

maximum  residual  bound  sup^g5;^^^.^  infmgVN  B{w,q;fj,)  decreases  monotonically  with  the  nested  enrichment  of 

Qn 

the  finite-element  space  and  the  hierarchical  enrichment  of  the  reduced-basis  space. 


5.  Numerical  Results 


5.1.  Case  Description 

We  consider  two  problems.  The  first  problem  is  a  parametrized  reaction-diffusion  problem  over  a  unit  square 
domain,  O  =  [0, 1]^.  The  governing  equation  is  given  by 

—V  •  (/iVrt)  -I-  u  =  1  in  n, 
u  =  0  on  dfl; 

here  jj,  G  T>  =  [0.01, 1]  is  the  reaction-coefficient.  The  solution  associated  with  the  two  extreme  values  of  the 
parameter  are  shown  in  Figure  1.  We  will  use  this  spatially  and  parametrically  smooth  problem  to  study 
the  convergence  behavior  of  our  minimum-residual  mixed  formulation  with  spatial  refinement  and  parametric 
refinement. 
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(d)'!i(Ai  =  l)  (e)pi(/i  =  l)  (f)p2(At  =  l) 


Figure  1.  Solution  to  the  reaction-diffusion  problem  for  two  different  parameter  values:  /r  = 

0.01  and  /i  =  1  (Poisson).  Note  that  u  is  the  primal  solution,  pi  is  the  first  component  of  the 
dual  solution,  and  p2  is  the  second  component  of  the  dual  solution. 

The  second  problem  is  a  parametrized  advection-diffusion  problem  over  a  L-shaped  domain,  =  [— 1, 1]^  \ 
[—1,0]^.  The  governing  equation  is  given  by 

-V  •  (Vu)  -f  V  •  {I5{p)u)  =  1  in  11, 
u  =  0  on  i9f2; 

the  advection  coefficient  is  given  by  /3{p)  =  (/i,  0)^  for  the  parameter  p  G  V  =  [0,  20].  The  solutions  associated 
with  the  two  extreme  values  of  the  parameter  are  shown  in  Figure  2.  Note  that,  due  to  the  presence  of  the 
reentrant  corner,  the  solutions  to  this  problem  belong  to  only  for  e  >  0. 

We  wish  to  construct  a  reduced-basis  approximation  that  achieves  the  residual  bound  tolerance  of  etoi  =  0.01 
for  any  p  G  T).  The  training  parameter  set  Strain  C  V  for  the  first  problem  consists  of  201  logarithmically 
equidistributed  points  between  [0.01,1];  the  training  parameter  set  for  the  second  problem  consists  of  201 
linearly  equidistributed  points  between  [0,20].  We  will  train  our  online  system  such  that  etoi  is  met  for  all 
p  G  Strain;  we  emphasize  that  the  residual  bound  with  respect  to  the  infinite-dimensional  function  space  may 
be  evaluated  for  any  p  that  is  not  necessarily  in  Strain,  but  the  residual  bound  tolerance  may  not  be  met  for 
M  ^  Strain-  We  employ  (7°  finite  elements  for  and  KT^  Raviart-Thomas  elements  for 

5.2.  Uniform  Spatio-Parameter  Refinement 

We  assess  the  behavior  of  the  minimum-residual  mixed  reduced-basis  method  under  uniform  refinement  in 
spatial  and  parameter  spaces.  Specifically,  we  employ  a  sequence  of  (nested)  spatially  uniform  finite- element 
meshes  and  a  sequence  of  (non-nested)  parametrically  equidistributed  snapshots. 
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(d)  =  20)  (e)pi(^  =  20)  (f)p2(/i  =  20) 


Figure  2.  Solution  to  the  L-shaped  advection-diffusion  problem  for  two  different  parameter 
values:  /i  =  0  (Poisson)  and  /i  =  20.  Note  that  u  is  the  primal  solution,  pi  is  the  first  component 
of  the  dual  solution,  and  p2  is  the  second  component  of  the  dual  solution. 

The  convergence  of  the  method  applied  to  the  reaction-diffusion  problem  with  the  finite-element  space  di¬ 
mension  M  is  shown  in  Figure  3(a);  the  convergence  of  the  method  with  the  reduced-basis  space  dimension  N 
is  shown  in  3(b).  In  Figure  3(a),  we  observe  that,  given  a  sufficient  reduced-basis  dimension  (i.e.  N  =  4),  the 
residual  bound  decays  at  the  rate  of  =  0{h?)  with  finite-element  refinement;  in  other  words,  we  achieve 
the  optimal  convergence  rate  for  P^-KT^  mixed  elements  in  norm  (which  is  equivalent  to  our  residual).  We 
also  note  that,  if  the  reduced-basis  space  is  not  sufficiently  rich,  then  the  convergence  of  the  residual  bound  with 
finite-element  refinement  stagnates.  In  Figure  3(b),  we  observe  that,  given  a  sufficient  finite-element  resolution 
(i.e.  Af  =  229377),  the  residual  bound  decays  exponentially  with  the  reduced-basis  refinement;  this  is  consistent 
with  the  theoretical  convergence  result  anticipated  for  the  smooth  parametric  manifold  [3].  We  also  note  that, 
if  the  finite-element  space  is  not  sufficiently  rich,  then  the  convergence  of  the  residual  bound  with  reduced-basis 
refinement  stagnates.  The  result  signifies  the  importance  of  improving  both  the  spatial  resolution,  through  the 
finite-element  enrichment,  and  the  parameter-space  resolution,  through  the  reduced-basis  enrichment.  In  the 
context  of  the  “standard”  reduced  basis  method  with  the  error  bounds  relative  to  the  finite-element  “truth” ,  this 
result  identifies  the  finite-element  fidelity  required  to  rigorously  justify  the  “truth”;  such  a  rigorous  validation 
is  in  general  non-existent  and  is  often  overlooked  in  the  “standard”  reduced-basis  procedure. 

We  now  assess  the  convergence  of  the  method  applied  to  the  advection-diffusion  problem  with  the  corner 
singularity.  Figure  4(a)  shows  the  convergence  of  the  residual  certificate  with  the  finite-element  space  dimension 
Af.  We  focus  the  case  with  a  sufficient  reduced-basis  resolution:  N'  =  6;  we  observe  that,  due  to  the  presence  of 
the  singularity,  the  asymptotic  convergence  rate  is  limited  to  despite  the  use  of  the  P^-KT^ 

elements.  Figure  4(b)  shows  that  we  indeed  need  N  =  172033  to  achieve  the  desired  residual  bound  of  10“^;  on  a 
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(a)  convergence  with  M  (b)  convergence  with  N 


Figure  3.  The  reaction-diffusion  problem.  Convergence  of  the  maximum  residual  bound  over 
the  training  set  Strain  C  T)  with  the  dimension  of  (a)  the  finite-element  space  and  (b)  the 
reduced-basis  space. 

smaller  finite-element  space,  we  cannot  achieve  the  residual  tolerance  regardless  of  the  reduced-basis  refinement 
due  to  the  lack  of  the  spatial  resolution.  We  note  that  the  required  number  of  the  finite-element  degrees  of 
freedom  is  quite  high  —  M  =  172033  —  due  to  the  presence  of  the  singularity  even  for  this  seemingly  simple 
problem.  If  we  wish  to  obtain  a  higher- fidelity  solution  that  meets  the  residual  bound  of,  say,  10“^,  then  the 
asymptotic  estimate  suggests  that  we  would  need  N  =  0(10®)  —  a  prohibitively  high  computational  cost. 

5.3.  Adaptive  Spatio-Parameter  Refinement 

We  now  invoke  the  spatio-parameter  Greedy  algorithm  to  solve  the  advection-diffusion  problem.  We  initialize 
the  algorithm  with  the  coarse  finite-element  mesh  shown  in  Figure  5(a)  that  consists  of  only  six  elements.  The 
first  parameter  selected  by  the  algorithm  is  /i  =  0.^  Since  the  finite-element  error  dominates  the  reduced-basis 
error  on  this  coarse  mesh,  the  algorithm  invokes  ten  steps  of  adaptive  mesh  refinement:  the  number  of  degrees  of 
freedom  increases  from  43  to  4495,  and  the  residual  bound  for  /i  =  0  reduces  from  0.24  to  0.0068.  We  emphasize 
that  our  formulation  provides  a  bound  even  on  the  unreasonable  coarse  initial  mesh,  since  the  bound  is  uniform 
and  not  asymptotic.  As  shown  in  Figure  6(a),  the  bound  decreases  at  the  rate  of  which  is  the  optimal  rate 
for  the  P^-MT^  discretization  despite  the  presence  of  the  singularity  thanks  to  the  adaptive  refinement.  (For  a 
theoretical  analysis,  see  for  instance  a  monograph  by  Schwab  [16].)  The  mesh  after  the  ten  steps  of  refinement, 
j-k=io ^  is  shown  in  Figure  5(b);  we  observe  a  localized  refinement  towards  the  singular  corner. 

Upon  meeting  the  residual  bound  tolerance  on  for  /i  =  0,  the  algorithm  continues  with  the  greedy 

sampling  of  the  parameter  space.  The  reduced-basis  sampling  of  the  parameter  space  identifies  that  the  residual 
bound  is  maximized  for  =  20.  The  algorithm  hence  solves  for  =  20)  and  (/i  =  20)  with  the  intention 

of  augmenting  Vat^i  =  span{M-^i®(^  =  0)}  and  Qn=i  =  span{p^i®(^  =  0)}.  However,  the  algorithm  identifies 
that  the  residual  bound  associated  with  ^  =  20  on  is  0.060  >  etoi-  The  algorithm  hence  invokes  two 

additional  steps  of  mesh  refinement  for  /i  =  20:  the  number  of  degrees  of  freedom  increases  from  4495  to  9521, 
and  the  residual  bound  for  =  20  reduces  from  0.060  to  0.0088.  The  error  decreases  rapidly  with  a  relatively 

^For  N  =  0,  the  residual  is  identical  for  all  S  X*;  the  algorithm  picks  the  first  occurrence,  =  0,  in  this  case. 
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M  =  ATy  +  Mq 


N 


(a)  convergence  with  M 


(b)  convergence  with  N 


Figure  4.  The  advection-diffusion  problem.  Convergence  of  the  maximum  residual  bound 
over  the  training  set  Strain  C  T)  with  the  dimension  of  (a)  the  finite-element  space  and  (b)  the 
reduced-basis  space. 


rf: 


(b)  rt 


(<:)  rt 


Figure  5.  finite-element  meshes:  (a)  initial  mesh;  (b)  adapted  mesh  after  refinement  at  /x  =  0; 
(c)  adapted  mesh  after  additional  refinement  at  /x  =  20. 


small  increase  in  the  number  of  degrees  of  freedom,  as  shown  in  Figure  6(b).  The  mesh  after  the  two  steps  of 
refinement,  is  shown  in  Figure  5(c);  we  observe,  compared  to  shown  in  Figure  5(b),  additional 

refinement  in  the  boundary  layer. 

Upon  meeting  the  residual  bound  tolerance  for  /x  =  20  —  and  also  for  /x  =  0  thanks  to  the  nested  finite- 
element  refinement  —  the  algorithm  continues  with  the  greedy  sampling  of  the  parameter  space.  This  time 
algorithm  completes  the  construction  of  the  reduced-basis  spaces  Vat^s  C  and  C  without 

requiring  additional  finite-element  enrichment.  The  greedy  convergence  over  the  parameter  domain  with  the 
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J\f  =  Afy  +  Mq 
(a)  /^  =  0  :  r^'==°  ^ 

Figure  6.  Convergence  histories 


for  the  adaptive  finite-element. 


(a)  parametric  variation  (b)  maximum  residual 


Figure  7.  Convergence  of  reduced-basis  approximation  on 


dimension  of  the  reduced-basis  space,  N,  is  shown  in  Figure  7.  With  the  adaptive  finite-element  snapshots 
each  of  which  meets  the  desired  residual  tolerance,  we  observe  the  exponential  convergence  with  N.  It  is  also 
worth  noting  the  considerable  reduction  in  the  number  of  finite-element  degrees  of  freedom,  Af,  achieved  by  the 
adaptive  finite-element  method  compared  to  the  case  with  uniform  meshes. 

6.  Summary  and  Perspectives 

We  propose  a  reduced  basis  method  for  parametrized  PDFs  certified  by  a  bound  on  the  dual  norm  of  the 
residual  computed  in  an  infinite-dimensional  functional  space.  The  residual  bound  builds  on  a  minimum-residual 
mixed  finite  element  method  and  the  associated  mixed  reduced  basis  approximation.  We  emphasize  that  our 
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method  constructs  reduced  basis  spaces  that  meet  the  residual  bound  tolerance  for  all  training  points  Strain  C  T), 
and  not  just  the  points  about  which  the  snapshots  are  collected.  In  addition,  in  the  online  stage,  we  may  assess 
if  the  reduced-basis  approximation  meets  the  required  tolerance  for  any  parameter  value  regardless  of  whether  it 
belongs  to  the  offline  training  set.  Hence,  the  greedy  algorithm  presented  here  is  fundamentally  different  from  a 
two-step  procedure  of  1)  invoking  the  standard  reduced-basis  method  that  computes  the  residual  relative  to  the 
finite-element  space  —  not  the  infinite-dimensional  function  space  —  to  identify  least-well-represented  parameter 
and  2)  invoking  an  adaptive  finite-element  method  in  computing  the  snapshots;  the  two-step  procedure  can  in 
general  only  provide  infinite-dimensional  residual  bound  only  for  the  snapshots. 

Because  the  proposed  minimum-residual  mixed  reduced  basis  method  yields  a  coercive  system  and  provides 
uniform  residual  bound  for  any  underlying  finite-element  discretization  and  any  reduced-basis  snapshots,  we 
may  explore  various  adaptive  strategies.  At  the  finite-element  level,  we  may  consider  more  general  hp  and 
anisotropic  adaptivity  in  the  physical  space  [16].  At  the  reduced-basis  level,  we  may  explore  hp  adaptivity  in 
the  parameter  space  as  considered  by  Eftang  et  al.  [6] ;  we  may  also  incorporate  finer-grain  mixing  of  snapshots 
computed  on  different  finite-element  spaces  as  considered  by  Steih  and  Urban  [17].  The  design  of  various  spatio- 
parameter  adaptive  schemes  —  enabled  by  the  proposed  offline-online  efficient  uniform  residual  bound  —  is  a 
key  to  effectively  treat  spatially  and  parametrically  complex  problems. 


Appendix  A.  Coercivity  of  the  Bilinear  Form 

We  prove  in  this  appendix  Proposition  2,  the  coercivity  of  the  bilinear  form  G(-,  •;  /r):  for  any  p  G  V, 

G{{v,q),{'v,q);h)  >  a{p){\\v\\l  +  Ugllg)  Vu  G  V,  Vg  e  Q, 


for 


l;(/i)  = 


Wp)y‘ 


(l  +  2||G(-; /x)||£(y_^2(Q))  -I-  2||A(- 


;  A^)ll£(v.(L2(n))rf; 


-1 


Proof.  For  notational  simplicity,  in  this  proof,  we  denote  the  inner  product  and  norm  over  O  without  explicit 
subscripts:  that  is  (•,•)  =  and  ||  •  ||  =  ||  •  ||L2(f2).  In  addition,  we  suppress  the  explicit  appearance  of 

the  parameter  p  in  various  forms  and  operators. 

By  way  of  preliminaries,  we  recall  that  our  bilinear  form  a(-,-)  is  inf-sup  stable  with  the  inf-sup  constant 
(3  >  0.  We  introduce  the  associated  supremizer:  S{w)  G  V  such  that 


=  a{w,v)  Vu  G  V; 


P  = 


inf  sup 


l{w,  v) 


l|ia||vlFllv 


inf 

wGV 


II^HIIv 

Ikllv 


note  that 
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We  then  note  that,  for  any  fc  €  M, 

G{{v,q),  {v,q))  =  ||C'(p)  +  V  •  gf  +  P(w,  Vp)  -  gf  +  ||g  •  n||^2(r„) 

=  ||C(p)f  +  2(C(n),  V  •  g)  +  II V  •  gf  +  \\A{v,  Vt;)f  -  2{Aiv,  Vp),  g)  +  ||gf 
+  Ik  •  ^Ili2(r„)  -  2A:(S'(n),  V  •  g)  -  2fc(VS'(v),  g)  +  2k{S{v),  q  ■  n)L2(r„) 

=  ||C(i;)f  +  2{C{v)  -  kS{v),\/-  g)  +  |1V  •  gf 
+  IIAf ,  Vp)f  -  2{A{v,Wv)  +  kWSiv),q)  +  ||gf 
+  Ik  •  ^Ili2(r„)  +  2{kS{v),  q  ■  n)L2(r„) 

=  ||V  •  g  +  (Civ)  -  kSiv))r  +  ||C(n )f  -  ||C(n)  -  )f 
+  ||g  -  {A{v,  Vv)  +  ))f  +  Pf ,  Vn)f  -  Pf,  Vv)  +  )f 

+  ||g  •  n  +  fc-5'(w)|||2(r„)  —  lk>5'('p)f  2(r„) 

=  ||V  •  g  +  (C(n)  -  fc5P)f  +  2k{C{v),  S{v))  -  f  )f 

+  ||g  -  Pf ,  Vp)  +  VfcS'(p))f  -  2k{A{v,Vv),V S{v))  -  fck|VS'(p)f 
+  ||g  •  n  +  fc-5'(w)||i2(r„)  —  fc^l|5'(n)||i2(r„) 

=  ||V  •  g  +  (C(n)  -  fc5P)f  +  ||g  -  P(n,  Vn)  +  fcV^(n))f  +  ||g  •  n  +  fc5(n)||i2(r„) 

+  2ka{v,S{v))  -  fc^||S'(i;)||v 
>  2ka{v,  S{v))  —  fc^||S'(p)||y 

=  2k\\S{v)\\l-k^S{v)rv 

=  k{2-k)\\Siv)\\l. 

Some  clarifications  are  in  order:  the  first  equality  follows  from  the  definition;  the  second  equality  follows  from 
the  Green’s  theorem;  the  fourth  equality  follows  from  completing  the  square;  the  sixth  equality  follows  from 
a{w,  v)  =  —{A{w,  w),  v)  +  {C{w),v)  Vw,  v  €  V  and  |k||y  =  ||  Vpf  +  Iff  +  lkllL2(r„)  €  V;  the  second  to  last 
equality  follows  from  a(w,v)  =  {S{w),v)\>  yw,v  G  V.  We  now  set  /c  =  1  to  obtain 

G{{v,q),{v,q))  >  Iff)!!?;- 
We  finally  appeal  to  the  definition  of  the  inf-sup  constant  to  obtain 

lkllv<;^lf(fllv<^G(Pg),Pg)). 

We  next  note  that 

II V  •  gf  <  2\\Civ)  +  V  •  gf  +  2||CPf  <  2G(P  g),  {v,  g))  +  2||Gf  fll?, 

<  2  (l  +  r^lPf  G((n,  g),  f ,  g))  Vp  e  V,  Vg  e  Q. 

We  similarly  note  that 

||gf  <  2PP  Vu)  -  gf  +  2PP  Vi;)f  <  2G(f ,  g),  {v,  g))  +  2pf  PIf 

<2(1+ P  ^PIl£(v,(L2(n))d)^  G((i;,  g),  (p,  g))  Vu  G  V,  Vg  e  Q. 
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The  combination  of  the  last  three  inequalities  yield 


Ikll?.  + 


la< 


4  +  /3  ^  (^1  +  2||C'||£(y_^2(n)) 


2|l^ll£(v,(L2(n))rf; 


J  G{{v,q),{v,q)) 
Vv  e  V,  yq  G  Q, 


which  is  the  desired  inequality. 
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